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■ Summary. At the end of their hves low mass stars such as our Sun lose most of 

their mass. The resulting planetary nebulae show a wide variety of shapes, from 
spherical to highly bipolar. According to the generalized interacting stellar winds 
model, these shapes are due to an interaction between a very fast tenuous outflow, 
and a denser environment left over from an earlier slow phase of mass loss. Previous 
analytical and numerical work shows that this mechanism can explain cylindrically 
symmetric nebulae very well. However, many circumstellar nebulae have a multipo- 
\ lar or point-symmetric shape. With two-dimensional calculations, Icke showed that 

fT^ ■ these seemingly enigmatic forms can be easily reproduced by a two-wind model 

in which the confining disk is warped, as is expected to occur in irradiated disks. 
Here, we present the extension to fully three-dimensional adaptive mesh refinement 
simulations of such an interaction. 

o 

1 Introduction 

O , 

|H , In the final phases of stellar evolution, low mass stars, such as our Sun, first 

C/3 i swell up and shed a dense, cool wind in the asymptotic giant branch (AGB) 

phase. This episode is followed by a fast, tenuous wind that is driven by the 
exposed stellar core, the future white dwarf. The planetary nebulae (PNe) 
\ resulting from this expulsion phase come in a wide variety of shapes, from 

j_j ■ spherical to highly bipolar. Some even have a multipolar or point-symmetric 

' shape. Balick [2] proposed that such forms arise due to an interaction between 

a slow disk-shaped inner AGB nebula and the fast 'last gasp' of the star. 
Analytical [12; 16] and numerical [39; 13; 28] work shows that this generalized 
interacting stellar winds (GISW) mechanism works very well (for an up-to- 
date review, see Balick & Frank [3]). Several scenarios for obtaining a disk 
around a PN exist, and it is in general assumed in the models that the shape 
of the dense gas around the star is a disk or a toroid. 

Icke [14] proposed that the point-symmetric shapes observed in a num- 
ber of PNe are formed in an interaction between a spherical stellar wind and 
a warped disk. It is prossible to produce such a warp around a single star. 



2 Erik-Jan Rijkhorst, Vincent Icke, and Garrelt Mellema 

through the combined effects of irradiation and cooUng [e.g. 34; 27] . Whereas 
Ickc's computations were restricted to a two-dimensional proof-of-principle, 
we now present a first series of fully three-dimensional hydrodynamic compu- 
tations of such a wind-disk interaction using the technique of adaptive mesh 
refinement (AMR). 

2 Point-symmetric nebulae 

Work on cylindrically symmetric nebulae showed [12; 13; 15] that sharply colli- 
mated bipolar flows are a frequent and natural by-product of the GISW model. 
However, many circumstellar nebulae have a multipolar or point-symmetric 
(i.e. antisymmetric) shape [38; 37]. The nebulae that are formed in the wind- 
disk interaction would naturally acquire the observed antisymmetry if the disk 
that confines the fast wind is warped, instead of symmetric under reflection 
about the equatorial plane. In binary stars such warps might be common, 
due to various torques in the system, but few point-symmetric nebulae have 
been proven to have binary cores. Several mechanisms have been proposed 
for warping a disk surrounding a single star. The most interesting one for our 
purposes invokes radiative instability [32; 17; 34; 27]. 

Livio & Pringle already proposed [23; 24] that the precession of warped 
disks might be responsible for point-symmetric nebulae. They proved con- 
clusively that the various physical scales for mass, accretion, luminosity and 
precession match the observations. The production of the nebulae proper they 
attributed to an unspecified 'jet' mechanism. Observations of many bipo- 
lar nebulae with 'ansae' (e.g. NGC 3242, NGC 7009) and 'FLIERS' (e.g. 
NGC 6751, NGC 7662) leave little doubt that jets are occasionally formed 
during the evolution of some aspherical PNe, probably in the late post-AGB 
phase, before the fast wind has switched on. But the nebulae presented by 
Sahai & Trauger [37] do not seem to resemble such shapes. While leaving open 
the possibility that jets may be responsible for additional structures, as in the 
case of the 'ansae', we show that the interaction between a warped disk and 
a spherically symmetric wind suffices. The lobes of point-symmetric nebulae 
[38; 37] look as if they were produced almost simultaneously. This is difficult 
in the case of a preccssing jet, which would make a corkscrew-like nebula of a 
type not readily apparent in post-AGB shells, although some objects do show 
features that are likely to be due to precession [38]. 

3 Radiation driven warping 

When an accretion disk is subject to external torques it may become unstable 
to warping [4; 32; 30] and when irradiated by a sufficiently luminous central 

star even an initially flat disk will warp [17; 34: 27: 26]. The difference in 
radiation pressure on slightly tilted annuli at different radii will induce the 
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warp. Essential is that the disk is optically thick for the stellar radiation and 

for its own cooling flux. The latter condition is the most restrictive, because 
a disk that is optically thin in the infrared will not suffice. This restricts the 
disks to a specific subclass of high density and low temperature. 

Analytical considerations lead to expressions for growth and precession 
rates and morphologies of the warp whereas numerical calculations including 
the effects of self-shadowing show that the non-linear evolution of the warp 
can produce highly distorted shapes, even with an inverted, counter-rotating 
inner disk region [35; 41]. Applications of warped disk theory range from 
active galactic nuclei [35; 26] to X-ray binaries [25; 41], protostellar disks [1], 
and PNe [23; 24]. Other mechanisms to produce warped disks, besides the 
radiatively driven one, are the wind driven [36] and the magnetically driven 
instability [20; 21]. 

Since we intend to study warped disks in PNe, we need a scenario to form 
accretion disks in these systems. Plausible scenarios include the coalescence 
of compact binaries [40], absorption of planetary systems, or the formation 
of a disk due to a main-sequence companion being out of equilibrium when 
emerging from a common envelope (CE) phase with a primary AGB star. In 
this case, the companion loses most of the mass it accreted during the CE 
phase which subsequently forms a disk around the primary [40] and which, in 
a later stage, can get radiatively warped when illuminated by the central star 
of the PN. 

For a PN the luminosity of the central star alone is sufficiently high to 
induce a radiation driven warp. Following Pringle [35], an expression for the 
radius Rcrit beyond which the disk is unstable to radiation driven warping 
is found from comparing the timescales of the viscous and radiation torques, 
leading to 

Rcrit = {2Tr/Af , (1) 

with the contant A defined by = l/4c~^G~^M~^Lj^~^tM~^. Here c is 
the speed of light, G the gravitational constant, rj = 1^2/ i^i is the ratio of the 
azimuthal to the radial viscosity, is the mass and L* the luminosity of 
the central star and Mace is the disk's accretion rate. We assumed a surface 
density = Macc/i^T^vi) [e.g. 33]. 

In a cartesian coordinate system, the warped disk surface is given by [34] 

cos (j) sin 7 -I- sin (j) cos 7 cos /3 \ 
cos (/) cos 7 -|- sin (/) sin 7 cos /3 , (2) 
— sin(^sin/3 J 

with local disk tilt angle /3(_R, 0), and orientation angle of the line of nodes 
7(i?, (/)). Here, R and (f) are the non-orthogonal radial and azimuthal coordi- 
nates respectively, pointing to the surface of the disk [cf. 34]. In our model 
calculations we adopt the case of a steady precessing disk with no growth 
and zero torque at the origin for which we have in the precessing frame that 
7 = A^fR and /3 = sin7/7, with the constant A defined by Eq. (1) [27]. 
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4 Time scales 



Since radiative cooling plays such an important role in our model (see the 
next section), we need to compare the cooling time scale to three other time 
scales related to the disk: the precession and growth time scales of the warp, 
and the shock passing time. The latter is defined as the time the expanding 
wind blown bubble takes to travel to the outer disk radius Rd- 

— 1/2 

Adopting a cooling function of the form A = AiTg and using the jump 
conditions for a strong shock, the cooling time of the shocked gas is given 
by [cf. 18; 19] tc = Cv^/pe, with Vg the shock speed, and pe the pre-shock 
environment density. When assuming a fully ionized gas of cosmic abundances, 
the constant C has a value of 6.0 x 10""^^ g cm~^ s^. 

Analytical relations for the radius Rait) and speed Vs{t) of the outer shock 
as functions of time can be derived [e.g. 22] and with these the shock passing 
time readily follows from sotting Rs{tsp) — Rd as 

^.p = (2/37rpe)'/'M-l/2«-l/2i^^ (3) 
The precession time scale of the disk is given by [27] 

tp = 48TT^cG^/^Ml^^L-^Rl/^Sd , (4) 

where we assumed Kcplcrian rotation. The time scale for the initial growth 
of the warp is of the same order. When we use Eq. (1) for the critical radius 
as a typical disk radius, we find that the different time scales scale as 

tp a MlL-^Mljf'Ed (5) 
tap oc pTM-'/W"MlL-^Mtcc^^ (6) 
t^ oc p-'/^MTvTM^'LlM-^lr' (7) 

We are quite limited in our choice of M*, L*, w^, and since values 

for these parameters arc strongly constrained by stellar evolution and wind 
models [e.g. 31; 7] but since the dependence of tp, i^p, and t^ on Alacc is so 
strong, a proper choice of this latter parameter leads to the desired proportion 
between the different time scales. 

For our calculations we used Myj = lO~*M0yr~^, Vyj = 200 kms~^, 
M* = 0.6Mq, = IO^Lq, Pe = lO-^^gcm-3, M„ec = IQ-^Mgyr-i, 
Ed = lgcm~^, and 77 = 1 resulting in Rcrit — 2AU, tp ~ 17yr, tsp ~ 0.4yr, 
tc — 10^*yr, and density contrast x = Pd/Pe — 300 with pd the disk density. 
So tc <C tsp <C tp showing that cooling will indeed be of importance and that 
we can safely ignore the disk's precession. 

5 Mechanism 

The mechanism behind the formation of the multipolar lobes s(h;ii in the 
simulations is as follows [see also 14]. As the central wind impinges on the 
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inner rim of the disk, a three-dimensional bow shock develops around it. This 

shock roughly has the shape of two 'cones' connected by there tips at the 
disk's center. The opening angle of the shock depends inversely on the Mach 
number of the wind. Imagining a two-dimensional cut through the developing 
bow shock, one sees that one branch flies off into space creating a lobe jutting 
out from the nebula, whereas the other slams into the concave side of the disk, 
scooping up disk material and thereby producing a set of smaller, unstable 
lobes (see Fig.l). Due to the cooling of the gas, the swept up shell is highly 
compressed and therefore thin, and the ram pressure from the wind directly 
drives the shell outwards, which are necessary ingredients for the bow shock 
to produce the lobes. Simulations of energy-driven flows never result in the 
pronounced point-symmetric lobes as seen in the calculations with cooling. 
When the density of the disk is not too high, the wind breaks through the 
concave part of the disk, producing another pair of lobes. 

Since the true bow shock is three-dimensional and the disk is warped, 
rotating the two-dimensional cut around the center of the disk shows that the 
concave side of the disk turns into the convex side and vice versa emphasizing 
the necessity of fully three-dimensional simulations of this interaction to truly 
understand the emerging point-symmetric structure. 

6 Numerical implementation 

We used the three-dimensional hydrocode Flash [11] to model the interaction 
between a spherical wind and a warped disk. This parallelized code implements 
block-structured AMR [6; 5] and a PPM type hydrosolver [42; 8]. 

Besides implementing the proper initial and boundary conditions, we also 
added our own cooling module in order to model radiative cooling using a 
cooling curve [10; 29]. This curve gives the energy loss rate as a function of 
temperature for a low density gas in collisional ionization equilibrium. The 
radiative losses were implemented through operator splitting and if the hydro 
timestep was larger than the cooling time, the former was subdivided into 
smaller steps when calculating the cooling. 

To construct the warped disk, Eq.(l) was combined with a constant 'wedge 
angle' 0^ and a proper value for A, i.e. was taken to be a few times 
Rcritj see Eq. (2). This disk was given a constant density which, through the 
density contrast x, resulted in a value for the environment density rig. The 
spherical wind was implemented as an inner boundary condition and given 
a 1/r^ density profile and a constant wind velocity Vyj. The pressure was 
calculated from an equation of state with a constant Poisson index 7. 

7 2D wind-disk simulations 

To check the implementation of the wind-disk interaction model into the AMR 
code described above, we repeated a number of the two-dimensional calcula- 



6 



Erik- Jan Rijkhorst, Vincent Icke, and Garrelt Mellema 











1 

■E 


KiiiSSS 

Ullllll 


■■■i;isllis9ill 

■■■^^■■■■{■■■■■■■■■iiBa 

■■■■■■■■■lilHHiHIIBi 

Mb:-' ^^9als!S!saiagHi5s 


■■■■ 


■ 

nil 
■lil 




■1 


■■■■■I 


■■■■■■■■■■■■■IIIBB 
■■■■■^H^IbBIIIIII 


nl 




1 




■i 


■ 


■ 


1 






■ 




■■ 




■ 










1 




1 


■ 


HI 




1 


1 






■! 




■1 ■■■■■ 






■ 






IIE 


■niiiiiiBI 


Uii!i 

H!!! 


iB 










IIB 


Ilia--! 


laiiSI 


■■■■■■■■■■I 


liiE^;: 

IIBBHI 
!!!■■ 


lie 
III 


■■■■1 

■■■■1 
■■■■ 


:■ 

M 

llil 

MUM 






■■ 
■■ 

■■ 

■ 


iHIH 
■■HSiiiJ 
BBIH 
■■■■■ 


IIHBI 
HBI 
HBI 

in 




aiai 




III 


■■■■IIIH 
■■■■■■III 
■■■■IIIII 

■ 


■ 




■■■■ 


mmmm 


■■ 




HI 

Hi 


■■■■ 
■■■■ 
■■■■ 
■■■■ 

sm 


Mill 
■■■■■■■I 

■■■■■■KE 

■mill 

■1 


III R^iilH 

mail 

ill 







Fig. 1. 2D example of a wind-disk interaction showing the evolving AMR grid 
structure superimposed on a plot of the logarithm of the density. Every square 
represents a grid of 8x8 cells. Five different levels of refinement are visible. The 
effective resolution of this simulation is 1024x1024 cells. 
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tions previously done by Icke [14] (see Fig.l). Because he used a different 
hydrosolvcr (FCT/LCD) and the outcome of the calculations strongly de- 
pend on turbulent processes in the gas, the simulations did not agree in every 
single detail but the overall point-symmetric morphologies were retrieved. All 
these simulations were nm using a small value for the Poisson index (7 = 1.1) 
resulting in 'momentum driven' bubbles. 

To see what happens when more realistic cooling is applied, we ran some 
simulations with the cooling curve module and a Poisson index 7 = 5/3. 
This showed that, apart from the production of the by now familiar point- 
symmetric lobes, the outer shell of swept up gas is thinner and unstable and 
developed into a number of smaller lobes merging with one another as the 
shell expanded. 

8 3D AMR Simulations 

Following the two-dimensional trial calculations, we ran wind-disk simulations 
in three dimensions on a cartesian grid with an effective resolution of up to 
512^ cells using five levels of refinement. Since we found in our two-dimensional 
calculations that simulations with cooling applied through a cooling curve do 
not result in a qualitatively different morphological outcome compared to 
calculations with a low value for the Poisson index (7 = 1.1), we opted not 
to use the cooling curve module for our three-dimensional simulations to save 
computational time. 

We used the following parameters: 7 = 1.1, rig = 5 x 10^ cm~^, x = 100) 
9d ~ 5°, and = 200 kms"-'^. In Fig. 2 we present a visualization of the 
three-dimensional shape of the swept up shell through isosurfaces at different 
viewing angles. Also shown are the corresponding synthesized Ha images, 
derived by projecting the three-dimensional data cube onto the plane of the 
sky. For this, we simply integrated the density squared along the line of sight 
and used this as a rough estimate for the emission. 

9 Conclusions 

Our computations show that the wind-disk interaction model in which the 
confining disk is warped results in a wide variety of point-symmetric shapes. 
Nebulae that show 'punched holes', such as NGC 7027 [9], are readily ac- 
commodated in this model. Other candidates for our model are probably 
NGC 6537 and, to a lesser extent, NGC 7026, in which the inner disk is 
still visible, while the outer nebula shows clear point-symmetric striictures. 
Also, the shapes of the (proto-)PNe He 2-155 and M 4-18 (and maybe 
He 2-47) can be explained with our model. As a further application, large- 
scale explosions in non-planar disks, such as might occur in active galax- 
ies, are expected to show similar patterns, provided that the disk mate- 
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Fig. 2. Isosurfaces of the density at the end of the wind-disk simulation as seen 
from different angles {left column) and the corresponding synthesized Ha images 
{right column). 
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rial can cool rapidly enough. Movies of these simulations can be found at 
http : //www . strw . leidenuniv . nl/AstroHydro3D/. 
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